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An inviscid two-dimensional fluid model with nonlinear dispersion that arises simultaneously in 
coarse-grained descriptions of the dynamics of the Euler equation and in the description of non- 
Newtonian fluids of second grade is considered. The scaling of the equilibrium states of this model 
for conserved energy and enstrophy retains the corresponding scaling for the Euler equations on the 
large scales and at the same time greatly deemphasizes the importance of small scales. This is the 
first clear demonstration of the beneficial effect of nonlinear dispersion in the model, and should 
highlight its utility as a subgrid model in more realistic situations. 



I. INTRODUCTION 

In 1998, Holm et al.,0 using the Euler-Poincare varia- 
tional formalism, proposed a model for the mean motion 
of ideal incompressible fluids. In this approach, the (re- 
duced) Lagrangian, which for the incompressible case is 
the kinetic energy, was modified from that for the Euler 
equation: 

l = \J |u| 2 dx (=E), 

to account for fluctuation energy of the velocity field in 
conjunction with the introduction of a fluctuation length 
scale an 

l = \J (|u| 2 + a 2 |Vu| 2 )dx (=£<*)• (1) 

The resulting "a-model" for the Euler equations is 

— + u • Vv - a 2 (Vu) T ■ V 2 u = -Vp 
ot 

V • u = 0, v = (1 - a 2 V 2 ) u (2) 

where, when a is set to zero, v = u, and the usual Euler 
equations are recovered. All other notation is standard. 
These equations are envisaged as modeling the flow of 
inviscid incompressible fluids at length scales larger than 
a. (For proof of existence and uniqueness of-solutions 
of (H), see Shkollem and Cioranescu & Girauhxl (viscous 
case).) p 

Rivlin and Ericksen,El in 1955, derived general consti- 
tutive laws of the differential type for an incompressible 
fluid, wherein at the first order, viscous Newtonian stress 
results (first grade fluids), while at the next order invis- 
cid, non-Newtonian, stress-strain relations appear (sec- 
ond grade fluids). Equations (||) are identically the equa- 
tions governing inviscid second-grade fluids, and where 
now a is a material property. Viscous and inviscid sec- 
ond grade fluid flows have since been studied from differ- 
ent viewpoints (e.g., see Dunn & FosdickO and references 
therein, and Cioranescu & Giraultcl). We also note that 
the variational formulation of (0) was already explicitly 
noted in Cioranescu & Ouazar.ET 

The new derivation of (^) has, however, renewed in- 
terest in them and besides spurring more mathematical 



work has stiHijftlated computational investigations of a- 
models( e '&-> EnI2P for the reason that the advection ve- 
locity, u, is obtained by a spatial-average of the advected 
field v (inversion of the Hclmholtz operator in (|J)). This 
results in a modification of the advective nonlinearity, 
the main nonlinearity of fluid dynamics, in such a way 
as to suppress mutual interactions between scales which 
are smaller than a (as can be seen, for example, in the 
untruncated version of (^]) below when \m\, \n\ > 2n/a). 
This modification is purely inviscid, and we will refer to 
it simply as nonlinear- dispersion in what follows. How- 
ever, with the exception of Nadiga & Shkoller,Ej compu- 
tational studies of a-models have always used additional 
viscous terms: For example, Chen et al.fl (1998) examine 
the applicability of a viscous a-model to model turbu- 
lent channel flow, and Chen et alfl (1999) explore the 
utility of a three-dimensional viscous a-model in provid- 
ing a subgrid model for fluid turbulence. While this is 
clearly the appropriate direction to pursue in the context 
of realistic applications, we think that studying purely 
inviscid a-models, although idealized, is important and 
will complement the studji of their viscous counterparts. 
In Nadiga and ShkollerjLlJ among other things, we pre- 
sented a series of two-dimensional numerical computa- 
tions comparing the solutions of Euler equations, Navier- 
Stokes equations, and an Euler-a model, and showed that 
the Euler-a model was able to reproduce the typical en- 
strophy decay characteristics of the Navier-Stokes equa- 
tions, but in a conservative setting. Presently, we address 
some statistical scaling aspects of the dynamics of such 
an Euler-a model to highlight its inviscid subgrid-scale 
modeling features. 

To better illustrate the effects of nonlinear- dispersion, 
the salient feature of all a-models, it suffices to consider 
in two-dimensions. In that case, it can be rewritten 
in the vorticity-streamfunction formulation as 

dq dq , t\ i l n 
d-t = d-t +J ^ q]= ° 

g=(l-a 2 V 2 )w, V 2 tP = lu, (3) 

where ip is the streamfunction, to is the vorticity, J is the 
Jacobian operator so that J[ip,q] = — + Tjjff^j an d 
again, when a is set to zero, q — lu, and the usual Euler 
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equations result. Equation (^) can be also be written as 



~dt 



a 2 V 2 ) 



J[j), (l - a 2 V 2 ) oj] = 0, V 2 ^ = uj, 



a form that highlights the modification to the J[?/>,w] 
nonlinear term of the Euler equations. Parenthetically, 
we note that in going to two-dimensions, wc lose analogs 
of three-dimensional processes like vorticity stretching, 
and therefore, fail to characterize the effect of nonlinear- 
dispersion on such processes. 

The kinetic energy E a (denoted by E when a = 0), as 
defined in ([!]), is an obvious constant of motion in both 
two and three dimensions. However, in two dimensions, 
unlike in three, the vorticity q (to when a — 0) of each 
fluid element is an inviscid constant (see (^|)), implying 
an infinity of conservation laws. In particular, enstrophy 
Z a , defined as 



Z a 



[(l-a 2 V 2 )w] 2 dx, 



(4) 



is a second conserved quadratic quantity. As before, 
when a = 0, we represent the conserved enstrophy by Z . 
(The domain integral of u • w or helicity is a quadratic 
quantity which is conserved in three dimensions, but 
which is identically zero in two dimensions.) 

The use of equilibrium statistical mechanical theo- 
ries (for (||) with a = 0) to better understand the 
inviscid dynamics of two-dimensional flows range from 
the two-constraint theory (see Kraichnan and Mont- 
gomerytLil and references therein) for finite truncations of 
the continuous system, to those based pn point vortices 
(again see Kraichnan and Montgomery^ and references 
thereiai-and their generalizations to continuous vorticity 
fieldsli3'llj which consider the infinity of conserved quanti- 
ties. In this article, we present the two-constraint theory 
for (|^) and verify the main results of the theory compu- 
tationally. Other than mentioning that there is already 
some numerical evidenceiS which seems to suggest that 
individual solutions of the Euler-a model (|J) may indeed 
follow predictions made for the behavior of the ensemble- 
averaged solutions of the Euler equations by the more 
complicated statistical theories we do not consider such 
theories any more in this short note. Also, since one 
may be tempted to point to the shortcomings of the two- 
constraint theory for the Euler equations before consid- 
ering the utility of such a theory for the Euler-a model, 
we wish to point out that the importance of this work lies 
primarily in the comparison of the results for the Euler-a 
model to the classical results for the Euler equations. In 
so doing, the effects of nonlinear- dispersion, and its bene- 
ficial numerical ramifications, arc clearly highlighted. At 
the risk of belaboring the point further, we reemphasize 
that in considering the simple two constraint theory, we 
are in no way suggesting that the behavior of the ensem- 
ble averaged solution of the a-model @ (or their slightly 
viscous counterparts) will follow this theory in more real- 
istic situations; the limitation of this theory in predicting 



large-scale coherent structures in the a — case is well 
knownjli and carry over to the nonzero a case. 

Furthermore, from a numerical point of view, invis- 
cid computations of (^) which conserve two quadratic 
invariants are fairly easily realizable and more common- 
place than the more involved multisimplcctic schemes 
which are required for conserving a larger number of 
constraints. Also, while state of the art schemes of the 
latter kind can handle only tens of modes (because of 
an N 3 log N scaling of computational work, where, N 
is the number of modest-a), there is no such restriction 
on schemes of the former kind. Examples of schemes 
which conserve just the energy and enstrophy invariants 
are Fourier-Galerkin truncations implemented as a fully 
dealiased pseudospectral spatial discreteization and the 
second-order finite difference spatial discreteization us- 
ing the Arakawa Jacobian. While we have done compu- 
tations with both these schemes and see no discrepancy 
between the results, we consider only the spectral dis- 
cretization in this article since the theory presents itself 
most naturally in this setting. 



II. TWO-CONSTRAINT STATISTICAL THEORY 
FOR (3) 

Let q x represent a discretization of q on a two- 
dimensional spatial grid, x, with 2A+1 equispaced points 
on each side. Let qu, where k is the set of all wave- 
vectors k — (k x ,k y ) denote the Fourier transform of 
q x . Although there are (2N + l) 2 A;-space grid points, 
since q x is real, not all of them are independent and 
q k = <J*L k - Therefore, there are only half as many fc-space 
grid points, and, k the set of all k is such that 

k = \k — {kxi ^y)? k m i n k X: ky ^ k max } • 

However, since each q k is a complex number, there are 
overall (2N + l) 2 degrees of freedom in q\. Consider the 
truncation of (ft) that is closed in q\: 



dt qk 



£ 

m-\-n=k 



m x n 



m\ 2 (1 + a 2 |m| 2 ) 



= 0. 



(5) 



Among the infinity of conservations for the continuous 
system (||) previously discussed, conservations (0) and 
([I|) are the only ones which survive for the truncated 
system (^|) , and may be expressed in terms of q k as 



E 



\Qk 



2^|£f(l + a 2 |*fr 



Z" 



lzZ^\ 2 - 

feek 



(G) 



This follows from the detailed conservation property of 
energy and enstrophy wherein each of these quantities is 
conserved in every triad interaction. 

Considering the dynamics of q\ under (|^), we work 
in the (2 A + l) 2 dimensional phase space. As a conse- 
quence of (S) satisfying a detailed Liouville theorem (see 



2 



Kraichnan and Montgomery^ and references therein), 
(||) also satisfies a Liouville theorem and the motion of 
$k in the truncated phase space is divergence freeJlil We 
can, therefore, define a stationary probability density, P, 
such that P rifcek ^Qk 1S the probability of finding the 
system within the ((2N + l) 2 dimensional) phase space 
volume Ilfcek^^ centered around (fk, and the ensemble 
average of any quantity O, a function of 5k, as 



(O) = [OPJI dq k 



(7) 



Next, a maximization of the information theoretic en- 
tropy s, defined in the usual fashion as 

S = -{bxP - 1) = - / (PlnP-P) J] dq k , 
J kek 

subject to constant ensemble-averaged energy and en- 
strophy, (E a ) and (Z a ) respectively, leads to 



P = aexp(-(3E a -jZ a ) 



(8) 



Here, ft (an inverse temperature associated with energy) 
and 7 (an inverse temperature associated with enstro- 
phy) are the Lagrange multipliers associated with the 
two constraints, and a is determined from 



J feek 



= f . 



Making use of (g) in (g) then leads to a factorization of 
the probability density: 



P = a Y[ CX P ( 



feek 



ft 



|fc| 2 (l + a 2 |£f) 



(9) 



The ensemble averaged two-dimensional spectral den- 
sity is then computed using (]7|) and (||) (after noting the 
expressions for the moments of a Gaussian) as 



1 



(U a (k)) = -< 



9k 



1 



1 



2'|fc| 2 (l + a 2 |fc |2 ^ 



4 ft + j\k\ 2 (1 + a 2 \k\ 2 ) 



Since (the isotropic) onc-dimensional spectra are more 
convenient for plotting, we define 

E a (\k\)= J2 ( U °U))> sothat ^ = ^^(|fc|). 

|fc|<|i|<|fc|+i |fc| 

In what follows, we drop the | • | sign on k and to avoid 
confusion, note that while E a represents the total con- 
served energy, E a {k), with a dependence on k, represents 
the corresponding one-dimensional spectrum. The one- 
dimensional spectrum E a {k) is then seen to scale with k 
as 



E a (k) 



k 



/3 + 7 fc 2 (l + a 2 fc 2 )' 



(10) 



with the above scaling being only approximate when the 
mode spacing is not small compared to k (as at small k). 

In (|l0|), since a is a given length scale, once the dis- 
cretization is fixed, expressions for the total energy and 
enstrophy of the given initial conditions provide two 
equations to solve for ft and 7. The equilibrium spectral 
scaling ( |To| ) is then seen to exhibit three regimes depend- 
ing on the values of the conserved energy and enstrophy 
as follows. If the minimum and maximum wavenumbers 
of the truncation are k m i n aivL-k max respectively, and if 
we define a mean wavenumbercJ of the initial conditions 



as, 




then, we can identify three regimes depending on the 
signs of ft and 7: 

• If the initial conditions are such that the mean 
wavenumber k\ is small: fe m ,„ < k\ < k a , then the 
temperature corresponding to energy is negative, 
while that corresponding to enstrophy is positive: 
-7*w(! + a2k Ln) < /3 < 0, 7 > 0; 



• If the mean wavenumber ki is medium: k a < k\ < 
fct,, then both temperatures are positive: ft > 0, 
7>0; 

• If the mean wavenumber k\ is large: k < k\ < 
k m ax, then the temperature corresponding to en- 
ergy is positive while the temperature correspond- 
ing to enstrophy is negative: ft > 0, —ft < 
lk 2 max {l + a 2 k 2 max ) <0. 

Here, k a , and fcf, are constants depending on the filter 
length a and the discretization: 



h 2 

■2 ^rnax 



k 2 



log 



k max (\ + a k min ) 

fcmin(l + & k max ) 



k 2 

■ 2 max 



V 2 k 4 
v min _j_ ^2 max 



Ul 1,2 , iL4 
max^min ' rnin 



(In the case of an infinite domain, the first of the above 
cases, ft < 0, cannot occur since k a =0.) 

Further, we can also compute the spectrum of the en- 
ergy conserved by the Euler equation (E) under the dy- 
namics of the Euler-a model. Noting that 



E = 



2 ^ 

kek 



\9k\ 



£f(l + a 2 |fc| 2 ) 2 ' 



(an extra factor (l+a 2 |fc| 2 ) in the denominator compared 
to the expression for E a ) and that E is not conserved 
for a / 0, the scaling of its one-dimensional spectrum, 
denoted simply by E(k), may be written ast£l 



E(k) 



(1 + a 2 k 2 ) {ft + jk 2 (1 + a 2 fc 2 )) ' 



(11) 
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III. DISCUSSION AND COMPUTATIONAL 
VERIFICATION OF RESULTS 

We devote the remainder of the article to a discussion 
of the scalings ( |l~0| ) and (|ll]) and their computational 
verification. First, when a is set to zero, in either ( |f(i| ) 
or (|Tl|), the classic result of Kraichnantil for the Euler 
equation: 



E(k) 



fc 



+ ~fk 2 



(12) 



is recovered, with the three regions corresponding to the 
different combination of signs for (3 and 7 now separated 
by values of the mean wavenumber k% corresponding to 
k a and kb, where k a and kb are given by 



kl 



log 



;.2 



k 2 4- k 2 

max ~ ^min 



As has been noted many times now0 for a = 0, there 
is no discontinuity of any sort in going from one region 
to the other among the three regions corresponding to 
different combinations of signs of [3 and 7. Therefore, for 
convenience, we first consider, in detail, the case (3 > 
and 7 > 0, and define 



fc? = 




lim fc° 



and note that 



K = K (1 



0(a 2 fc 2 )) . 

Furthermore, fc* can be shown0 to be of the order of k\. 
(Thus, for simplicity in what follows, one may use fci, 
fc" , fc* , and fc" interchangably, or represent all of them 
by ki.) For the Euler solutions, we have from ( |lG| ) with 
a = 0, the large scales and small scales (with respect to 
fc*) behaving asymptotically as 



E(k) ~ k, 
E{k) ~ k-\ 



'max i 



(13) 



implying equipartition of E at large scales and equipar- 
tition of Z at small scales. (When J: < i„ 7/c 2 <C (3 in 
(|l2|) and when k > ky_jk 2 > [3 in (Q.) When a is not 
zero, however, from (J10|), one easily sees the analogous 
E a - and Z Q -equipartition results to be respectively 



E a (k) ~ k, 
E a (k) ~ fc" 3 , 



fc* ^ fc !^ k max . 



(14) 



This implies that nonlinear- dispersion in (^) acts in such 
a way as to preserve the Euler scaling of dynamics 



at the large scales while at the same time greatly 
deemphasizing the importance of small scales. 

Asymptotic scalings arising from ( fil| ) , for the noncon- 
served energy for a ^ 0: 



rain — ^ ^- fv * : 



E(k) - k, k 
E(k) - fc" 5 , k* < fc < fc 



(15) 



further reinforce this xesult. Finally, we note that for 
a = 0, it is well knownEj (and easy to see from (|To|) ) that 
when k max — » 00, energy diverges logarithmically and 
enstrophy diverges quadratically. However, when a^0, 
one can see from (|l^) that when k max — > 00, energy is not 
divergent and that the enstrophy Z a is quadratically di- 
vergent. Nevertheless, in Nadiga and Shkollcr,EJ we show 

that it is the dynamics of the non-conserved enstrophies 

i 2 



<ix which are actu- 
the 



i/^xand I^(l-a 2 V 2 ) 

ally interesting. While the former does not diverge 
latter diverges only weakly (logarithmically). 

We have carried out a series of computational exper- 
iments on a doubly periodic two-dimensional domain, 
wherein an ensemble of initial conditions were evolved 
under @ for different values of a until statistical equili- 
bration was achieved. Initial conditions, similar to those 
in Fox and Orsag,E3 were obtained by choosing ampli- 
tudes for wavenumbers in the band 50 < |fc| < 51 (zero 
elsewhere) from a zero-mean normal distribution of ran- 
dom numbers. The variance was scaled in such a wajz. 
that for the different values of a, the conserved energytj 
had the same value: E a = E. The mean wavenumber k\ 
for this set of initial conditions corresponds to about 50.5, 
and for the resolution chosen below, 28.5 < k a < 43.5, 
and kb > 60.1. With this setup, besides letting both j3 
and 7 be positive, we can realize both the energy and 
enstrophy equipartition regimes in the same experiment. 

A fully dealiased pseudospectral spatial discretization 
was used with k m i n — 1, and k max — N = 85, and 
the nonlinear terms were computed in physical space us- 
ing 256 grid points in each direction. A nominally fifth- 
order, adaptive time stepping, Runge-Kutta Cash-Karp 
algorithmnJ was used for time integration. Energy was 
conserved to better than 1 in 10 5 and enstrophy to better 
than 1 in 10 4 over the entire duration of the runs. 

In Fig. [l], we plot the (instantaneous) spectrum E a (fc) 
against fc, on a log-log scale for four different values of 
a: 0, 0.05, 0.10, and 0.15, corresponding to between 
to 2.4 percent of the domain size. The spectrum for the 
Euler case (a = 0) is offset by a decade so as to not 
clutter the figure, and slopes of 1, —1, and —3 are drawn 
for reference. The scalings, ( |l3| ) for a — 0, and ( p"4| ) for 
nonzero a, are clearly verified at both large and small 
scales, in Fig. [I], with in fact a cleaner (but identical) 
scaling of the large scales for nonzero a. Furthermore, 
the three spectra corresponding to the nonzero-a cases 
seem to collapse onto a single curve in Fig. [l[ This col- 
lapse may be explained by noting that almost all the 
energy is contained in the low-fc modes (denoted below 
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by the set of wavenumbers kg) and almost all enstrophy 
is contained in the high-/c modes (denoted below by the 
set of wavenumbers kz). This leads to a leading order 
expression for (3 which is independent of a: 

and one for 7 which is inversely proportional to a : 
7 (a)~]Tfc/(a 2 fc^). 

kz 

This in turn implies that the spectra ([n]) should be al- 
most independent of a, except for a small intermediate 
range of k. 
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FIG. 1. Plot of E a (k) vs. k for q=0 (offset by a decade), 
and a = 0.05, 0.10, & 0.15. Low-fc scaling is identical for zero 
or nonzero a, but high-fc scaling is much steeper for a / 0. 
The theoretical asymptotic slopes are drawn for reference. 

The above collapse of the nonzero-a spectra onto a sin- 
gle curve seems to suggest that the actual value of a is 
not very important as long as its value is in a certain 
range. This, however, is not true as can be seen from the 
structure of the E(k) spectrum plotted in Fig. || for the 
same four cases. The theoretical scaling for this spec- 
trum is given in ([ll]) and the asymptotic scalings are 
given in (|l5|). Although the small scale behavior is as 
expected, it is clear that with increasing a, the struc- 
ture of the large scales is being substantially modified. 
(Scaling (|l5|) at large scales, can obviously be realized 
by increasing the number of modes, but that is not our 
intent here; we are examining the effect of a at fixed 
resolution.) Therefore, a small value of a is indicated. 
In such a case, the nonlinear- dispersion of the a-model 
is highly beneficial at small-scales while the large-scale 
distortion is minimal. Considering that the minimally 
resolved length scale in these computations corresponds 
to about 0.074, one may conclude that a should be of 



that order. That is to say, besides their use in describ- 
ing second-grade fluids, a-models in general and (|J) in 
particular should be useful as a subgrid model in under- 
resolved computations. 
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FIG. 2. Plot of Eik) vs. k for a = 0, 0.05, 0.10, & 0.15. 
Spectrum falls off much faster (k~ 5 ) for a 7^ compared 
to fc _1 for a = 0, but the large-scales are also increasingly 
changed with increasing a. 

These conclusions are also borne out in numerical ex- 
periments corresponding to the fluid-dynamically more 
interesting case wherein the temperature associated with 
energy is negative (/3 < 0). As mentioned earlier, such 
a case is obtained when the initial conditions are chosen 
with energy and enstrophy such that ki < k a . (As be- 
fore, 28.5 < k a < 43.5 for the different values of a for 
the discretization chosen.) In such a case, there is a cop-, 
densation of energy on to the low modes of the systemtLil 
resulting in large scale structures (necessarily coherent). 
However, the enstrophy equipartition scaling of the spec- 
tra discussed previously are unchanged. This gives us an 
opportunity to better test the extent of distortion of the 
low wavenumber (coherent) modes due to increasing a. 
While various aspects of the negative temperature f$se 
for nonzero a are considered in Nadiga and Shkoller,EHl in 
the spirit of this article, we presently consider only the 
spectral distortion to the structure of the low wavenum- 
ber (coherent) modes. In Fig. ||, where we plot the spec- 
trum E{k) versus k again for a — 0, 0.05, 0.10, & 0.15, 
and now where k\ w 2. For this case, only one realization 
(for each value of a) is considered and the spectrum cor- 
responds to a long time average, for good measure, taken 
after the system has reached statistical equilibrium. For 
this suite of runs, energy was conserved to machine preci- 
sion while enstrophy was conserved to about 0.3% for the 
entire duration of the runs considered. The steep slope 
of —5, for the small scales, for nonzero a (compared to a 
slope of —1 for a — 0) is again verified and more impor- 
tantly, for the case of a = 0.05, the low-mode structure 
(up to k = 10) is almost identical to the case of a = 0. 
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This is clearly not the case for the two other values of 
the filter length, a, which are greater than the smallest 
resolved scale of the computation. 
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FIG. 3. Negative temperature case. Plot of E(k) vs. k for 
a = 0, 0.05, 0.10, & 0.15. Small-scale spectrum falls off much 
faster (k~ 5 ) for a / compared to fc" 1 for a = 0. While the 
large scale distortion is minimal for a=0.05, it is appreciable 
for the other two cases with larger a. 
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